home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / haeberli / libgutil / stream.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  16KB  |  806 lines
  1. /*
  2.  * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17. /*
  18.  *    stream - 
  19.  *        Convert an image into lines.
  20.  *
  21.  *                Paul Haeberli - 1987
  22.  *
  23.  */
  24. #include "image.h"
  25. #include "obj.h"
  26. #include "vect.h"
  27. #include "stream.h"
  28.  
  29. static statecoords();
  30. static nextstate();
  31. static statesame();
  32. static unmarkcorners();
  33. static marksharpcorners();
  34. static badpoop();
  35. static lineintersect();
  36.  
  37. object *imgtoshape();
  38.  
  39. #define BOXTOL        (1.5*stoutscale)    /* for ratherstraight */
  40. #define CORNER        (1.2)             /* for marksharpcorners */
  41. #define SHORTSTEP    (3.0*stepsize)         /* for fixcorners */
  42. #define SHARP        (1.0)
  43. #define NOTSHARP     (0.0)
  44. #define SKIPPOINTS    (1)
  45.  
  46. #define XBRIDGE        (1<<8)
  47. #define YBRIDGE        (2<<8)
  48.  
  49. #define getpixel(x,y)        ((*(stfb+((y)*stxsize)+(x)))&0xff)
  50. #define setpixel(x,y,val)    *(stfb+((y)*stxsize)+(x)) = (val)
  51. #define setpixelbit(x,y,val)    *(stfb+((y)*stxsize)+(x)) |= (val)
  52. #define clrpixelbit(x,y,val)    *(stfb+((y)*stxsize)+(x)) &= ~(val)
  53. #define getpixelbit(x,y,val)    (*(stfb+((y)*stxsize)+(x)) & (val))
  54.  
  55. static short *stfb;
  56. static int stxsize, stysize;
  57. static float stoutscale;
  58.  
  59. object *fbtoshape();
  60. static float toslope();
  61. static float fakeatan();
  62. static float pntdist();
  63. static float pntslope();
  64. static float slopedelta();
  65. static point *pntinterp();
  66. static object *followpath();
  67. static object *resampleobj();
  68. static object *smoothobj();
  69. static object *filterobj();
  70. static object *fixcorners();
  71.  
  72. #define DX    0    /* X axis */
  73. #define DY    1    /* Y axis */
  74. #define MINUS    0    /* minus direction */
  75. #define PLUS    1    /* plus direction */
  76.  
  77. typedef struct pathstate {
  78.     int x;        /* origin x and y */
  79.     int y;
  80.     int axis;        /* DX or DY */
  81.     int direction;    /* PLUS or MINUS */
  82. } pathstate;
  83.  
  84. pathstate statetab[4][3] = {
  85.     {   { -1, 0, DY, MINUS },    /* from DX MINUS */
  86.     { 0, -1, DX, MINUS },
  87.     { 0, 0, DY, PLUS }     },
  88.  
  89.     {    { 0, 1, DY, PLUS },    /* from DX PLUS */
  90.     { 0, 1, DX, PLUS },
  91.     { -1, 1, DY, MINUS },    },
  92.  
  93.     {    { 0, 0, DX, PLUS },    /* from DY MINUS */
  94.     { -1, 0, DY, MINUS },
  95.     { 0, -1, DX, MINUS },    },
  96.  
  97.     {    { 1, -1, DX, MINUS },    /* from DY PLUS */
  98.     { 1, 0, DY, PLUS },
  99.     { 1, 0, DX, PLUS },       },
  100. };
  101.  
  102. #define FIX_SHARP    1
  103. #define FIX_ISLONG    2
  104. #define FIX_DELETE    4
  105.  
  106. typedef struct workpnt {
  107.     float x, y;
  108.     float slope;
  109.     short flags;
  110. } workpnt;
  111.  
  112. object *imgtoshape(name,thresh,func,stepsize)
  113. char *name;
  114. int thresh, func;
  115. float stepsize;
  116. {
  117.     int xsize, ysize;
  118.     short *fb;
  119.     object *obj;
  120.  
  121.     sizeofimage(name,&xsize,&ysize);
  122.     fb = (short *)shortimagedata(name),
  123.     obj = fbtoshape(fb,xsize,ysize,thresh,func,stepsize);
  124.     free(fb);
  125.     return obj;
  126. }
  127.  
  128. static zaphigh()
  129. {
  130.     int x, y;
  131.     short *sptr;
  132.     int val;
  133.  
  134.     sptr = stfb;
  135.     for(y=0; y<stysize; y++) 
  136.     for(x=0; x<stxsize; x++) 
  137.         *sptr++ &= 0xff;
  138.     val = getpixel(0,0);
  139.     for(x=0; x<stxsize; x++) 
  140.     setpixel(x,0,val);
  141.     for(x=0; x<stxsize; x++) 
  142.     setpixel(x,stysize-1,val);
  143.     for(y=0; y<stysize; y++) 
  144.     setpixel(0,y,val);
  145.     for(y=0; y<stysize; y++) 
  146.     setpixel(stxsize-1,y,val);
  147. }
  148.  
  149. object *fbtoshape(myfb,xsize,ysize,thresh,func,stepsize)
  150. short *myfb;
  151. int xsize, ysize, thresh, func;
  152. float stepsize;
  153. {
  154.     int x, y;
  155.  
  156.     stfb = myfb;
  157.     stxsize = xsize;
  158.     stysize = ysize;
  159.     if(stxsize>stysize)
  160.     stoutscale = 1.0/stxsize;
  161.     else
  162.     stoutscale = 1.0/stysize;
  163.     zaphigh();
  164.     for(y=1; y<stysize; y++) {
  165.     for(x=1; x<stxsize; x++) {
  166.         if(getpixel(x,y) < thresh) {
  167.         if(getpixel(x-1,y) >= thresh) 
  168.             setpixelbit(x,y,XBRIDGE);
  169.         if(getpixel(x,y-1) >= thresh) 
  170.             setpixelbit(x,y,YBRIDGE);
  171.         } else {
  172.         if(getpixel(x-1,y) < thresh) 
  173.             setpixelbit(x,y,XBRIDGE);
  174.         if(getpixel(x,y-1) < thresh) 
  175.             setpixelbit(x,y,YBRIDGE);
  176.         }
  177.     }
  178.     }
  179.     return followpath(thresh,func,stepsize);
  180. }
  181.  
  182. static object *followpath(thresh,func,stepsize)
  183. int thresh, func;
  184. float stepsize;
  185. {
  186.     int x, y;
  187.     pathstate initstate, curstate;
  188.     float xcoord, ycoord;
  189.     int count;
  190.     object *objlist = 0;
  191.     object *path = 0;
  192.  
  193.     count = 0;
  194.     for(y=1; y<stysize; y++) {
  195.     for(x=1; x<stxsize; x++) {
  196.         if(getpixelbit(x,y,XBRIDGE)) {
  197.         if(getpixel(x,y) >= thresh) {  /* XXX change */
  198.             initstate.x = x;
  199.             initstate.y = y;
  200.             initstate.axis = DX;
  201.             initstate.direction = MINUS;
  202.         } else {
  203.             initstate.x = x;
  204.             initstate.y = y;
  205.             initstate.axis = DX;
  206.             initstate.direction = PLUS;
  207.         }
  208.         curstate = initstate;
  209.         statecoords(&curstate,thresh,&xcoord,&ycoord);
  210.         path = newobj();
  211.         addpoint(path,newpnt(stoutscale*xcoord,stoutscale*ycoord,0.0));
  212.         while(1) {
  213.             nextstate(&curstate);
  214.             if(curstate.axis == DX) 
  215.             clrpixelbit(curstate.x,curstate.y,XBRIDGE);
  216.             else
  217.             clrpixelbit(curstate.x,curstate.y,YBRIDGE);
  218.             if(++count == SKIPPOINTS) {
  219.             statecoords(&curstate,thresh,&xcoord,&ycoord);
  220.             addpoint(path,newpnt(stoutscale*xcoord,stoutscale*ycoord,0.0));
  221.             count =0;
  222.             }
  223.             if(statesame(&curstate,&initstate)) 
  224.             break;
  225.         }
  226.         clrpixelbit(initstate.x,initstate.y,XBRIDGE);
  227.         statecoords(&initstate,thresh,&xcoord,&ycoord);
  228.         path = filterobj(path,func,stepsize);
  229.         if(path) {
  230.             path->type = OBJ_LOOP;
  231.             if(!objlist) 
  232.             objlist = path;
  233.             else 
  234.             addobject(objlist,path);
  235.         }
  236.         }
  237.     }
  238.     }
  239.     return objlist;
  240. }
  241.  
  242. static object *filterobj(obj,func,stepsize)
  243. object *obj;
  244. int func;
  245. float stepsize;
  246. {
  247.     object *path;
  248.     int n;
  249.  
  250.     if(func&RESAMPLEFUNC) {
  251.     path = resampleobj(obj,stepsize);
  252.     freeobj(obj);
  253.     if(!path)
  254.        return 0;
  255.     obj = path;
  256.     }
  257.  
  258.     n = unmarkcorners(obj);
  259.     if(n<3) {
  260.     freeobj(obj);
  261.     return 0;
  262.     }
  263.     marksharpcorners(obj);
  264.  
  265.     if(func&SMOOTHFUNC) {
  266.     path = smoothobj(obj);
  267.     freeobj(obj);
  268.     obj = path;
  269.     }
  270.  
  271.     if(func&FIXCORNERFUNC) {
  272.     path = fixcorners(obj,stepsize);
  273.     freeobj(obj);
  274.     obj = path;
  275.     }
  276.  
  277.     n = unmarkcorners(obj);
  278.     if(n<3) {
  279.     freeobj(obj);
  280.     return 0;
  281.     }
  282.     return obj;
  283. }
  284.  
  285. static statesame(s1,s2)
  286. pathstate *s1, *s2;
  287. {
  288.     if(s1->x != s2->x)
  289.     return 0;
  290.     if(s1->y != s2->y)
  291.     return 0;
  292.     if(s1->axis != s2->axis)
  293.     return 0;
  294.     if(s1->direction != s2->direction)
  295.     return 0;
  296.     return 1;
  297. }
  298.  
  299. static nextstate(s)
  300. pathstate *s;
  301. {
  302.     pathstate *next;
  303.     int x, y;
  304.     int major, minor;
  305.  
  306.     major = (s->axis<<1)+s->direction;
  307.     x = s->x;
  308.     y = s->y;
  309.     for(minor=0; minor<3; minor++) {
  310.     next = &statetab[major][minor];
  311.     if(next->axis == DX) {
  312.         if(getpixelbit(x+next->x,y+next->y,XBRIDGE)) {
  313.         s->x += next->x;
  314.         s->y += next->y;
  315.         s->direction = next->direction;
  316.         s->axis = next->axis;
  317.         return;
  318.         }
  319.     } else {
  320.         if(getpixelbit(x+next->x,y+next->y,YBRIDGE)) {
  321.         s->x += next->x;
  322.         s->y += next->y;
  323.         s->direction = next->direction;
  324.         s->axis = next->axis;
  325.         return;
  326.         }
  327.     }
  328.     }
  329.     badpoop(0);
  330. }
  331.  
  332. static statecoords(s,thresh,rx,ry)
  333. pathstate *s;
  334. int thresh;
  335. float *rx, *ry;
  336. {
  337.     int pixval0, pixval1;
  338.     float p0, delta;
  339.     int x, y;
  340.  
  341.     x = s->x;
  342.     y = s->y;
  343.     pixval0 = getpixel(x,y);
  344.     if(pixval0 < thresh) {
  345.     if(s->axis == DX) {
  346.         pixval1 = getpixel(x-1,y);
  347.         if(pixval1 >= thresh) {
  348.         *rx = x + 0.5+(pixval0-thresh)/((float)(pixval1-pixval0));
  349.         *ry = y + 0.5;
  350.         } else 
  351.         badpoop(1);
  352.     } else {
  353.         pixval1 = getpixel(x,y-1);
  354.         if(pixval1 >= thresh) {
  355.         *rx = x + 0.5;
  356.         *ry = y + 0.5+(pixval0-thresh)/((float)(pixval1-pixval0));
  357.         } else 
  358.         badpoop(1);
  359.     }
  360.     } else {
  361.     if(s->axis == DX) {
  362.         pixval1 = getpixel(x-1,y);
  363.         if(pixval1 < thresh) {
  364.         *rx = x + 0.5+(pixval0-thresh)/((float)(pixval1-pixval0));
  365.         *ry = y + 0.5;
  366.         } else 
  367.         badpoop(1);
  368.     } else {
  369.         pixval1 = getpixel(x,y-1);
  370.         if(pixval1 < thresh) {
  371.         *rx = x + 0.5;
  372.         *ry = y + 0.5+(pixval0-thresh)/((float)(pixval1-pixval0));
  373.         } else 
  374.         badpoop(1);
  375.     }
  376.     }
  377. }
  378.  
  379. static badpoop(n)
  380. int n;
  381. {
  382.     fprintf(stderr,"imgtoshape: bad poop %d\n",n);
  383.     exit(1);
  384. }
  385.  
  386. static object *fixcorners(obj,stepsize)
  387. object *obj;
  388. float stepsize;
  389. {
  390.     workpnt *shape;
  391.     workpnt *p, *np, *start, *end;
  392.     point *pnt;
  393.     int val, i, j, m, npoints;
  394.     object *robj;
  395.     float delta, dx, dy, z;
  396.     vect p1, p2, p3, p4, pinter;
  397.  
  398. /* count the points */
  399.     npoints = 0;
  400.     pnt = obj->points;
  401.     while(pnt) {
  402.     npoints++;
  403.     pnt = pnt->next;
  404.     }
  405.  
  406. /* put them into a handy data structure */
  407.     shape = (workpnt *)mymalloc(npoints*sizeof(workpnt));
  408.     pnt = obj->points;
  409.     p = shape;
  410.     for(i=0; i<npoints; i++) {
  411.     p->x = pnt->x;
  412.     p->y = pnt->y;
  413.     p->flags = 0;
  414.     if(pnt->z == SHARP) {
  415.         p->flags |= FIX_SHARP;
  416.     }
  417.     pnt = pnt->next;
  418.     p++;
  419.     }
  420.  
  421. /* add length and slope info */
  422.     for(i=0; i<npoints; i++) {
  423.     p = shape+i;
  424.     np = shape+((i+1)%npoints);
  425.     dx = p->x - np->x;
  426.     dy = p->y - np->y;
  427.     if(sqrt(dx*dx+dy*dy) > SHORTSTEP) {
  428.         p->flags |= FIX_ISLONG;
  429.     }
  430.     p->slope = toslope(np->x-p->x,np->y-p->y);
  431.     }
  432.  
  433. /*
  434.  * find place where we have long edge, then some number of sharp points
  435.  * followed by another long edge. 
  436.  */
  437.     for(i=0; i<npoints; i++) {
  438.     p = shape+i;
  439.     if(p->flags & FIX_ISLONG) {
  440.         j = 0;
  441.         while(1) { 
  442.             p = shape+((i+j+1)%npoints);
  443.         if((p->flags & FIX_SHARP) == 0)
  444.            break;
  445.         if(p->flags & FIX_ISLONG)
  446.            break;
  447.         j++;
  448.         }
  449.         if(j>0 && j<3 && (p->flags & FIX_SHARP) && 
  450.                            (p->flags & FIX_ISLONG)) {
  451.         start = shape+((i+0)%npoints);
  452.         end = p;
  453.         delta = slopedelta(start->slope,end->slope);
  454.         if((delta>0.5) && (delta<3.5)) {
  455.             for(m=1; m<=j; m++) {
  456.             p = shape+((i+m+1)%npoints);
  457.             p->flags |= FIX_DELETE;
  458.             }
  459.             p = shape+((i+0)%npoints);
  460.             p1.x = p->x;
  461.             p1.y = p->y;
  462.             p = shape+((i+1)%npoints);
  463.             p2.x = p->x;
  464.             p2.y = p->y;
  465.             p = shape+((i+j+1)%npoints);
  466.             p3.x = p->x;
  467.             p3.y = p->y;
  468.             p = shape+((i+j+2)%npoints);
  469.             p4.x = p->x;
  470.             p4.y = p->y;
  471.             lineintersect(&p1,&p2,&p3,&p4,&pinter);
  472.             p = shape+((i+1)%npoints);
  473.             p->x = pinter.x;
  474.             p->y = pinter.y;
  475.         }
  476.         }
  477.     }
  478.     }
  479.  
  480. /* make a linked list object to return */
  481.     robj = newobj();
  482.     p = shape;
  483.     pnt = 0;
  484.     for(i=0; i<npoints; i++) {
  485.     if(!(p->flags & FIX_DELETE)) {
  486.         if(p->flags & FIX_SHARP) 
  487.         z = SHARP;
  488.         else 
  489.         z = NOTSHARP;
  490.         if(!pnt) {
  491.         robj->points = pnt = newpnt(p->x,p->y,z);
  492.         } else {
  493.             pnt->next = newpnt(p->x,p->y,z);
  494.         pnt = pnt->next;
  495.         }
  496.     }
  497.     p++;
  498.     }
  499.     free(shape);
  500.     return robj;
  501. }
  502.  
  503. static point *ratherstraight(start,firstpoint,n)
  504. point *start, *firstpoint;
  505. int n;
  506. {
  507.     int count;
  508.     point *end, *p;
  509.     float dist;
  510.     float a, b, c;     /* for  ax+by+c == 0 */
  511.     float mag;
  512.  
  513.     end = start;
  514.     for(count=0; count<n; count++) {
  515.     if(end->next) {
  516.         end = end->next;
  517.         if(end->z == SHARP && (count != n-1)) 
  518.         return 0;
  519.     } else {
  520.         if(count == n-1)
  521.         end = firstpoint;
  522.         else
  523.         return 0;
  524.     }
  525.     }
  526.     if(n == 1)
  527.     return end;
  528.     else {
  529.        a = end->y - start->y;
  530.        b = end->x - start->x;
  531.        mag = sqrt(a*a+b*b);
  532.        a /= mag;
  533.        b /= mag;
  534.        c = (a*start->x)-(b*start->y);
  535.        p = start->next;
  536.        while(p!=end && p) {
  537.        dist = (a*p->x)-(b*p->y)-c;
  538.        if(dist<0.0)
  539.            dist = -dist;
  540.        if(dist>BOXTOL)
  541.            return 0;
  542.        p = p->next;
  543.        }
  544.        return end;
  545.     }
  546. }
  547.  
  548. static object *smoothobj(obj)
  549. object *obj;
  550. {
  551.     object *cobj;
  552.     point *pnt, *cpnt, *lastpnt, *npnt;
  553.     float curslope;
  554.     int n;
  555.  
  556.     cobj = newobj();
  557.     pnt = obj->points;
  558.     cpnt = newpnt(pnt->x,pnt->y,pnt->z);
  559.     cobj->points = cpnt;
  560.     while(pnt) {
  561.     n = 1;
  562.     lastpnt = 0;
  563.     while(1) {
  564.         if(npnt = ratherstraight(pnt,obj->points,n)) {
  565.         if(npnt==obj->points) 
  566.             return cobj;
  567.         if(npnt->z == SHARP) {
  568.             cpnt->next = newpnt(npnt->x,npnt->y,npnt->z);
  569.             cpnt = cpnt->next;
  570.             pnt = npnt;
  571.             break;
  572.         }
  573.         lastpnt = npnt;        
  574.         n++;
  575.         } else {
  576.         if(lastpnt) {
  577.             cpnt->next = newpnt(lastpnt->x,
  578.                         lastpnt->y,lastpnt->z);
  579.             cpnt = cpnt->next;
  580.             pnt = lastpnt;
  581.             break;
  582.         } else {
  583.             pnt = pnt->next;
  584.             break;
  585.         }
  586.         }
  587.     }
  588.     }
  589.     return cobj;
  590. }
  591.  
  592. static object *resampleobj(obj,stepsize)
  593. object *obj;
  594. float stepsize;
  595. {
  596.     object *cobj;
  597.     point *pnt, *cpnt;
  598.     float curslope;
  599.     float totallength;
  600.     float step;
  601.     float need, used;
  602.  
  603.     stepsize *= stoutscale; 
  604.     cobj = NULL;
  605.     totallength = 0;
  606.     pnt = obj->points;
  607.     while(pnt) {
  608.     if(pnt->next)
  609.         pnt->z = pntdist(pnt,pnt->next);
  610.     else
  611.         pnt->z = pntdist(pnt,obj->points);
  612.     totallength += pnt->z;
  613.     pnt = pnt->next;
  614.     }
  615.     if(totallength>(2*stepsize)) {
  616.     step = totallength/(int)((totallength/stepsize));
  617.     cobj = newobj();
  618.     pnt = obj->points;
  619.     cpnt = newpnt(pnt->x,pnt->y,0.0);
  620.     cobj->points = cpnt;
  621.     need = step;
  622.     used = 0.0;
  623.     while(pnt) {
  624.         if(((1.0-used)*pnt->z) < need) {
  625.         need -= (1.0-used)*pnt->z;
  626.         used = 0.0;
  627.         pnt = pnt->next;
  628.         } else {
  629.         if(pnt->next)
  630.             cpnt->next=pntinterp(pnt,pnt->next,used+need/pnt->z);
  631.         else 
  632.             cpnt->next=pntinterp(pnt,obj->points,used+need/pnt->z);
  633.         cpnt = cpnt->next;
  634.         used += need/pnt->z;
  635.         need = step;
  636.         }
  637.     }
  638.     if(pntdist(cobj->points,cpnt)<(step/10.0))
  639.         cobj->points = cobj->points->next;
  640.     }
  641.     return cobj;
  642. } 
  643.  
  644. static unmarkcorners(obj)
  645. object *obj;
  646. {
  647.     point *pnt;
  648.     int n;
  649.  
  650.     pnt = obj->points;
  651.     n = 0;
  652.     while(pnt) {
  653.     pnt->z = NOTSHARP;
  654.     pnt = pnt->next;
  655.     n++;
  656.     }
  657.     return n;
  658. }
  659.  
  660. static marksharpcorners(obj)
  661. object *obj;
  662. {
  663.     point *pnt, *npnt, *nnpnt, *nnnpnt;
  664.     float slope1, slope2;
  665.  
  666.     pnt = obj->points;
  667.     if(pnt->next)
  668.     npnt = pnt->next;
  669.     else
  670.     npnt = obj->points;
  671.     if(npnt->next)
  672.     nnpnt = npnt->next;
  673.     else
  674.     nnpnt = obj->points;
  675.     if(nnpnt->next)
  676.     nnnpnt = nnpnt->next;
  677.     else
  678.     nnnpnt = obj->points;
  679.     while(pnt) {
  680.     if(slopedelta(pntslope(pnt,npnt),pntslope(nnpnt,nnnpnt))>CORNER) {
  681.         npnt->z = SHARP;
  682.         nnpnt->z = SHARP;
  683.     }
  684.     pnt = pnt->next;
  685.     npnt = nnpnt;
  686.     nnpnt = nnnpnt;
  687.     if(nnpnt->next)
  688.         nnnpnt = nnpnt->next;
  689.     else
  690.         nnnpnt = obj->points;
  691.     }
  692. } 
  693.  
  694. static float slopedelta(s1,s2)
  695. float s1, s2;
  696. {
  697.     float delta;
  698.  
  699.     if(s1<s2) 
  700.     delta = s2-s1;
  701.     else
  702.     delta = s1-s2;
  703.     if(delta>4.0)
  704.     delta = 8.0-delta;
  705.     return delta;
  706. }
  707.  
  708. static float pntslope(p1,p2)
  709. point *p1, *p2;
  710. {
  711.     if(p2) 
  712.     return toslope(p2->x-p1->x,p2->y-p1->y);
  713.     else
  714.     return 100.0;
  715. }
  716.  
  717. static float pntdist(p1,p2)
  718. point *p1, *p2;
  719. {
  720.     float dx, dy;
  721.  
  722.     dx = p1->x - p2->x;
  723.     dy = p1->y - p2->y;
  724.     return sqrt(dx*dx+dy*dy);
  725. }
  726.  
  727. static point *pntinterp(p1,p2,param)
  728. point *p1, *p2;
  729. float param;
  730. {
  731.     float x, y, z;
  732.     float a, b;
  733.  
  734.     a = param;
  735.     b = 1.0-a;
  736.     x = b*p1->x + a*p2->x;
  737.     y = b*p1->y + a*p2->y;
  738.     return newpnt(x,y,0.0);
  739. }
  740.  
  741. static float fakeatan(a,b) 
  742. float a, b;
  743. {
  744.     a = a/b;
  745.     if(a>0.5773502) 
  746.     return 0.666666 + (0.333333/(1.0-0.5773502))*(a-0.5773502);
  747.     else
  748.     return a*(0.666666/0.5773502);
  749. }
  750.  
  751. static float toslope(dx,dy)
  752. float dx, dy;
  753. {
  754.     int flipx, flipy;
  755.     float val;
  756.  
  757.     if(dx<0) {
  758.     flipx = 1;
  759.     dx = -dx;
  760.     } else
  761.     flipx = 0;
  762.     if(dy<0) {
  763.     flipy = 1;
  764.     dy = -dy;
  765.     } else
  766.     flipy = 0;
  767.     if(dy<dx)
  768.     val = fakeatan(dy,dx);
  769.     else
  770.     val = 2.0-fakeatan(dx,dy);
  771.     if(flipx)
  772.     val = 4.0-val;
  773.     if(flipy)
  774.     val = 8.0-val;
  775.     return val;
  776. }
  777.  
  778. static lineintersect(p1,p2,p3,p4,pinter)
  779. vect *p1, *p2, *p3, *p4, *pinter;
  780. {
  781.     float dx, dy, mag;
  782.     float x3, y3;
  783.     float x4, y4;
  784.     float xc, yc;
  785.     float dx2, dy2;
  786.  
  787.     dx = p2->x - p1->x;
  788.     dy = p2->y - p1->y;
  789.     mag = sqrt(dx*dx+dy*dy);
  790.     dx /= mag;
  791.     dy /= mag;
  792.  
  793.     p3->x -= p1->x;
  794.     p3->y -= p1->y;
  795.     p4->x -= p1->x;
  796.     p4->y -= p1->y;
  797.     x3 = dx*p3->x + dy*p3->y;
  798.     y3 = dx*p3->y - dy*p3->x;
  799.     x4 = dx*p4->x + dy*p4->y;
  800.     y4 = dx*p4->y - dy*p4->x;
  801.  
  802.     xc = x3 + ((y3*(x4-x3))/(y3-y4));
  803.     pinter->x = p1->x + dx*xc;
  804.     pinter->y = p1->y + dy*xc;
  805. }
  806.